Effects of Disorder on Synchronization of Discrete Phase-Coupled Oscillators 
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We study synchronization in populations of phase-coupled stochastic three-state oscillators char- 
acterized by a distribution of transition rates. We present results on an exactly solvable dimer as 
well as a systematic characterization of globally connected arrays of M types of oscillators {M = 2, 
3, 4) by exploring the linear stability of the nonsynchronous fixed point. We also provide results for 
globally coupled arrays where the transition rate of each unit is drawn from a uniform distribution of 
finite width. Even in the presence of transition rate disorder, numerical and analytical results point 
to a single phase transition to macroscopic synchrony at a critical value of the coupling strength. 
Numerical simulations make possible the further characterization of the synchronized arrays. 
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I. INTRODUCTION 

Synchronization in populations of phase-coupled non- 
linear stochastic oscillators, and the corresponding emer- 
gence of macroscopic coherence, appear pervasively in a 
tremendous range of physical, chemical, and biological 
systems. As a result, the general subject continues to 
be studied intensely in applied mathematics and theo- 
retical physicsfll, 12, S 0, Q. Since the pioneering work 
of Kuramoto [31] , emergent cooperation in these systems 
has been investigated from a myriad of perspectives en- 
compassing both globally and locally coupled, stochastic 
and deterministic, and large and small systems. And 
while Kuramoto's canonical model of nonlinear oscilla- 
tors, whose use is widespread because of its close kin- 
ship to the normal form describing general phase oscilla- 
tions, has proven spectacularly successful for character- 
izing synchronization, simple, phenomcnological models 
of synchronization have also proven useful in a variety 
of new contexts 0, 0, d, Q , most notably the character- 
ization of emergent synchronization as a nonequilibrium 
phase transition. 

We have shown [S, 9] that a model of three-state iden- 
tical phase-coupled stochastic oscillators is ideally suited 
for studying the nonequilibrium phase transition to syn- 
chrony in locally-coupled systems, owing in large part to 
its numerical simplicity. The utility of these studies rests 
on the well-established notion of universality, that is, on 
the contention that microscopic details do not determine 
the universal properties associated with the breaking of 
time-translational symmetry that leads to a macroscopic 
phase transition. Statistical mechanics is thus enriched 
by simplistic, phenomcnological models (the Ising model 
being the most ubiquitous example) whose microscopic 
specifics are known to be, at best, substantial simplifi- 
cations of the underlying quantum mechanical nature of 
matter, but whose critical behavior captures that of more 
complex real systems. In this spirit, our simple tractable 



model captures the principal features of the syncrhoniza- 
tion of phase-coupled oscillators. In the globally coupled 
(mean field) case our model undergoes a supercritical 
Hopf bifurcation. With nearest neighbor coupling, we 
have shown that the array undergoes a continuous phase 
transition to macroscopic synchronization marked by sig- 
natures of the XY universality class (Tol . [Tij . including 
the appropriate classical exponents (3 and v, and lower 
and upper critical dimensions 2 and 4 respectively. 

In this paper we focus on globally coupled arrays and 
expand our earlier studies to the arena of transition rate 
disorder. We start with a slightly modified version of 
our original model (explained below), in which identical 
synchronized units are governed by the same transition 
rates as individual uncoupled units. Then, in the spirit 
of the original Kuramoto problem ^] , we explore the oc- 
currence of synchronization when there is more than one 
transition rate and perhaps even a distribution of tran- 
sition rates among the phase-coupled oscillators. In par- 
ticular, we explore the conditions (if any) that lead to a 
synchronization transition in the face of a transition rate 
distribution, discuss the relation between the frequency 
of oscillation of the synchronized array and the transition 
rates of individual units, and explore whether or not the 
existence of units of different transition rates in the cou- 
pled array may lead to more than one phase transition. 

Our paper is organized as follows. In Sec. |TT] we de- 
scribe the model, including the modification (and associ- 
ated rationale) of our earlier scenario even for an array 
of identical units. In Sec. IIIII we present results for a 
dimer composed of two units of different intrinsic tran- 
sition rates. While there is of course no phase transi- 
tion in this system, it is instructive to note that there 
is a probability of sychronization of the two units that 
increases with increasing coupling strength. Section IIVI 
introduces disorder of a particular kind, useful for a num- 
ber of reasons that include some analytical tractability. 
Here our oscillators can have only one of J\f distinct tran- 
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sition rates, where is a small number. We pay partic- 
ular attention to the dichotomous case, Af = 2. This 
simple disordered system reveals some important general 
signatures of synchronization. We also consider the cases 
Af — 3 and Af ~ 4, but find that the Af — 2 case already 
exhibits most of the interesting qualitative consequences 
of a distribution of transition rates. In particular, we are 
able to infer the important roles of the mean and variance 
of the distribution. In Sec. |V] we generalize further to a 
uniform finite-width distribution of transition rates and 
explore this inference in more detail. Section IVII sum- 
marizes our results and poses some questions for further 
study. 

II. THE MODEL 

Our point of departure is a stochastic three-state model 
governed by transition rates g (see Fig. [T]), where each 
state may be interpreted as a discrete phase [1, Q. 
The unidirectional, probabilistic nature of the transitions 
among states assures a qualitative analogy between this 
three-state discrete phase model and a noisy phase oscil- 
lator. The linear evolution equation of a single oscillator 
is dP{t)/dt = MP{t), where the components Pi{t) of the 
column vector P{t) = {Pi{t) P2{t) P^it))'^ {T denotes 
the transpose) are the probabilities of being in states 
I = 1, 2, 3 at time t, and 
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The system reaches a steady state for P^ — Pj* — P3 = 
1/3. The transitions i i + 1 occur with a rough pe- 
riodicity determined by g; that is, the time evolution 
of our simple model qualitatively resembles that of the 
discretized phase of a generic noisy oscillator with the 
intrinsic eigenfrequency set by the value of g. 

To study coupled arrays of these oscillators, we couple 
individual units by allowing the transition rates of each 
unit to depend on the states of the units to which it is 
connected. Specifically, for N identical units we choose 
the transition rate of a unit ly from state i to state j as 




FIG. 1: Tliree-state unit witli transition rates g. 




FIG. 2: (Color online) Phase space surface of the steady state 
probability VX that both units of a dimer are in the same 
state, indicating perfect synchronization, for 72 = 0.5 and a 
range of 71 and a. 
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FIG. 3: The steady state probability Va that both units of 
a dimer are in the same state for a range of fi and a. Top: 
contour image. Bottom: a single curve for a = 2.2. In the 
latter case, it is clear that as /i rises, synchronization rapidly 
decreases. 



where S is the Kronecker delta, a is the coupling param- 
eter, g is the transition rate parameter, n is the number 
of oscillators to which unit f is coupled, and Nk is the 
number of units among the n that are in state k. Each 
unit may thus transition to the state ahead or remain in 
its current state depending on the states of its nearest 



3 



neighbors. In our earlier work we considered the globally 
coupled system, n = N — 1, and also nearest neighbor 
coupling in square, cubic, or hypercubic arrays, n — 2d 
{d — dimensionality). Here we focus on the globally cou- 
pled array. 

Our previous work used a slightly different form of the 
coupling (12), with Ni^i Ni. While the differences in 
these details do not in any way affect the characteriza- 
tion of the synchronization transitions, that earlier choice 
was numerically advantageous because it led to a phase 
transition at a lower critical value Oc of the coupling con- 
stant (cc = 1.5 in the globally coupled array) than other 
choices. A lower coupling in turn facilitates numerical 
integration of equations of motion because the time step 
that one needs to use near the phase transition must be 
sufficiently small, dt <C e~'^/g. However, that earlier 
coupling choice brought with it a result that is undesir- 
able in our present context (but was of no consequence 
before). In our earlier model, as the units become in- 
creasingly synchronized above the transition point, the 
average transition rate of a cluster becomes substantially 
dependent on the value of a; specifically, the transitions 
and cluster oscillation frequency slow as a is increased 
due to an exponential decrease in the transition proba- 
bility. To cite an explicit example, consider a small sub- 
system composed of units which are all in the same state 
at time t (that is, a cluster of units which are perfectly 
synchronized). The previous form of the coupling yields 
an exponentially small transition rate in this case, and 
hence the oscillation frequency of this microscopic clus- 
ter approaches zero for high values of a. Since here we 
specifically wish to analyze the effects of transition rate 
disorder, it is desirable to deal with a model in which 
the average transition rate of identical synchronized units 
depends only on their intrinsic transition rate parameter 
and not on coupling strength. The form ([2]) reduces sim- 
ply to the constant g when the coupled units are perfectly 
synchronized. While the critical coupling in this new ver- 
sion is higher than in our earlier model and hence is nu- 
merically less efficient, no other features of the synchro- 
nization transition are affected. This in fact supports 
the desired insensitivity of the interesting macroscopic 
features of the model to microscopic modifications. 

For a population of iV — > oo identical units in the mean 
field (globally coupled) version of this model we can re- 
place Nk/N with the probability Pk, thereby arriving 
at a nonlinear equation for the mean field probability, 
dP{t)/dt = M[p[t)]P{t), with 
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FIG. 4: The steady state probability Va that both units of 
a dimer are in the same state for (71,72) = (0.5, 1.5) as a is 
increased (solid line) . The points represent simulation results, 
where Va is measured as the fraction of time that both units 
are fully synchronized. 



that cross the imaginary axis at Oc = 3, indicative of a 
Hopf bifurcation at this value. Note that the oscillation 
frequency of the array at the critical point as given by 
the imaginary parts of the eigenvalues is = 735/2. A 
more detailed analysis [1, 0, [l2| shows the bifurcation to 
be supercritical. 
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Normalization allows us to eliminate P^it) and obtain a 
closed set of equations for Pi (t) and P2 (t) . We can then 
linearize about the fixed point (PijPj*) = (1/3,1/3), 
yielding a set of complex conjugate eigenvalues which de- 
termine the stability of this disordered state. Specifically, 
we find that 2\±/ g = (a — 3) ± i^/?>{\ -t- a), eigenvalues 



FIG. 5: The two pairs of complex conjugate eigenvalues for 
the dichotomously disordered system, M = 1. Top panel: 
(A-,A1). Bottom panel: (A+,AI5_). The coupling constant is 
in the range 2.9 < a < 4.5, and the transition rate parameters 
are chosen to be (71,72) = (1,3). In the bottom panel the 
critical value of a is Oc « 3.95. 



4 



III. DIMER 



Consider first the simplest "disordered array," namely, 
a mutually coupled dimer where one unit is characterized 
by g = 7i and the other by 5 = 72. In terms of the states 
(phases) Si and 5*2 of units 1 and 2, there are 9 possible 
dimer states, {Si, S2) = (1, 1), (1, 2), . . . , (3, 3), but it is 
not necessary to seek the ensemble distributions for all 
of these states in order to decide whether or not the two 
units are synchronized. We can directly write an exact re- 
duced linear evolution equation for the 3 states A, B, and 
C, where A corresponds to any situation where both units 
are in the same state [that is, (5*1, 5*2) = (1, 1), (2, 2), and 
(3,3)], state B corresponds to a situation where unit 1 
is one state "ahead" [(5*1, 5*2) = (2,1), (3, 2), and (1,3)], 
and state C corresponds to a situation where unit 2 is 
one state "ahead" [{Si,S2) = (1,2), (2, 3), and (3,1)]. 
The evolution equation for these states is the closed lin- 
ear set 



dr{t)/dt:^AV{t), 



(4) 



with V{t) the time dependent probability column vector 
(VAit) Vsit) Vcit)f and 

-71-72 &72 bji \ 

A^ \ 71 -&"^7i - ^72 b-^l2 , (5) 

72 b~^ji -671-6^^72/ 

and where we have introduced the abbreviation 

6 = e^ (6) 

This evolution equation is easy to derive from the defi- 
nition of the coupling, Eq. For example, when the 
system is in state A, it can either go to state B, which 
happens when unit 1 jumps ahead with transition rate 
71, or it can go to state C, which happens when unit 2 
jumps ahead with transition rate 72. Similarly, when the 
system is in state B, it can either jump to state A (when 
the lagging unit transitions forward) with transition rate 
&72 or jump to state C (when the leading unit transitions 
forward) with transition rate &~^7i. 

With normalization, Eq. (jH) becomes a 2-dimensional 
equation having steady state solution 



Klf + &'7i72 + 7I) 



K2^,2 



^ 7i + 72(71 + 72) 



(7) 



(1 + 6 + &2)(^2 + ^2) ^ (2 -f 63)7172 ■ 



The eigenvalues of the two-dimensional matrix obtained 
from A after implementing normalization have negative 
real parts for all positive values of the parameters a, 71, 
and 72, indicating that the fixed points given by Eq. ([7]) 
are stable. Hence, the system asymptotically tends to 
this steady state solution. We are particularly inter- 
ested in V^, the probablity for the system to be synchro- 
nized. In terms of the single relative width parameter 
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FIG. 6: (Color online) Upper panel: Stability boundary for 
the dichotomously disordered system. The contour ReA+ — 
is plotted in (71,72,(1) space. This contour indicates the crit- 
ical point, where the Hopf bifurcation occurs and the disor- 
dered solution becomes unstable. The region above the con- 
tour represents the synchronized phase. Lower panel: Stabil- 
ity boundary in terms of relative width parameter. 



width/mean, 

2(71 -72) 
(71 +72) 

(—2 < /i < 2), this probability is 



(8) 



^ 4(2 + b^) 



(2 + 6) 



1 



^ 6(2 + 26-6^) 
4(2-f 62)(2 + 6)/ 



(9) 



The probability of synchronization for a dimer of identi- 
cal units {^JL = 0) is thus V\ = 6/(2 -f 6) = e"/(2 + e''). 
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which increases with increasing couphng. This is the 
maximal syncrhonization; it is easy to ascertain that 
decreases with increasing as one would anticipate. 
The full behavior of V'X as a function of the various pa- 
rameters is shown in Figs. [2]|4l The gradual increase 
in synchronization probability with increasing coupling 
turns into a sharp transition as a function of a in the 
infinite systems to be considered below. The decreased 
synchronization probability when the frequencies of the 
two units become more dissimilar (increasing /i) will also 
be reflected in the dependence of the critical coupling on 
transition rate parameter disorder. 




IV. Af DIFFERENT TRANSITION RATES 

Next we consider globally coupled arrays of oscilla- 
tors that can have one of J\f different transition rate 
parameters, g = ju, u — l,...,Af. To arrive at a 
closed set of mean field equations for the probabilities 
we again go to the limit of an infinite number of oscil- 
lators, N oo. However, we must do so while pre- 
serving a finite density of each of the Af types of oscil- 
lators. The probability vector is now 3A/'-dimensional, 

Pit) = (Pim ^2,71 -P3,7i •■■ Pl,'f^^ P'2,iu P^,iuV ■ The 
added subscript on the components of P{t) keeps track of 
the transition rate parameter. Explicitly, the component 
f'i,7„ is the probability that a unit with transition rate 
parameter 5 = 7^ is in state i. The mean field evolution 
for the probability vector is the set of coupled nonlinear 
differential equations dP{t)/dt = M^r[P{t)]P{t), with 



/M 



M^[P{t)] = 



71 

M 







... 
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(10) 



Here 



-512 (7«) 9zi{lu) 
512 (7«) -523 (7«) 

ff23(7«) -53i(7«), 



(11) 



and 



7«exp 



L fc=l 



P, 



(12) 

The function fi^k) is the fraction of units which have a 
transition rate parameter g = ^k- 

As before, probability normalization allows us to re- 
duce this to a system of 2J\f coupled ordinary differen- 
tial equations. It is interesting to compare this setup 
with that of the original Kuramoto problem with noise, 
where a continuous frequency distribution is introduced 
and the governing equation is a nonlinear partial differen- 
tial equation [the Fokker-Planck equation for the density 
p{0,uj,t)] 13]. The discretization of phase in our model 




FIG. 7: Upper panel: Critical coupling ac as a function of fi 
for a dichotomous (A/" = 2) array of globally coupled oscilla- 
tors. The two curves represent the exact relationship (lower 
curve) and the small fi approximation (upper curve), respec- 
tively. Lower panel: The frequency of synchronous oscillation 
at the transition point. Lower curve is the approximation as 
predicted by ImA+, upper curve is the exact result. 



results instead in a set of 2J\f coupled nonlinear ordinary 
differential equations. 

While it is nevertheless still difficult to solve these 
equations even for small Af, we can linearize about the 
disordered state P{t) ~ (1/3 1/3 . . . 1/3)^ and arrive at 
a 2Af X 2J\f Jacobian of the block matrix form 



/Ji(7i) ^2(71) ^2(71) 

^2(72) Jl(72) J2(72) 



J 



J2(7l)\ 

^2(72) ^ 



VJ2(7Jv) J2{jn) ... J2{7n) Ji{in)J 
The blocks Ji{g) and J2{g) are given by: 

Mg) = 



-2g -g - ag/N\ 
-ag/N) -g + ag/Nj 



and 



J2{9) = 







-ag/N 
ag/N ag/N 



(13) 
(14) 

(15) 



While we explore this in more detail below only for small 
Af, we note that in general the Jacobian p3|) has Af pairs 
of complex conjugate eigenvalues, only one pair of which 
seems to have a real part that becomes positive with in- 
creasing coupling constant a. This implies that there is 
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a single transition to synchrony even in the presence of 
the transition rate disorder that we have introduced here. 
We go on to confirm this behavior for TV = 2, 3, and 4. 



A. Two transition rate parameters 

For the JV = 2 case, the four eigenvalues 
(A+, A^, A_, Al) of the Jacobian can be determined ana- 
lytically. We find 



^^'^^ = i [a — 6 ± i?(a, ^) cos (C(a, ^))] , 



71+72 8 
ImA± _ 1 

71+72 8 

where 



%/3(a + 2) ± B{a, ^) sin (C(a, ^)) 



(16) 



B{a, fi) = V2 [a^ - 6aV^ + 3/(a^ + 3^ 



,1/4 



C(a,/i) 



■ tan 



-V3(a^ - (q + 3)aj^) ' 
a2 + 3(a - l)/i2 



(17) 



Aside from an overall factor (71 +72), Eqs. ^6]) depend 
only on the relative width variable as defined in Eq. ([5]), 
and therefore the critical coupling Oc depends only on /i. 
As illustrated in Fig. [51 one pair of eigenvalues crosses the 
imaginary axis at a critical value a = Oc, but the other 
pair shows no qualitative change as a is varied. While 
this figure shows only the particular transition rate pa- 
rameter values (71,72) — (1,3), the qualitative features 
of these eigenvalues remain similar for the entire range of 
positive parameters. The upper panel of Fig. [S] depicts 
the contour ReA-|- = in (71 , 72 , a) space; this contour 
represents the critical surface and thus separates the syn- 
chronous and disordered phases. 

The critical coupling is the value of a at which RcA+ = 
(RcA_ does not vanish for any a). It is easy to ascer- 
tain that ImA+ does not vanish at Uc, so that the critical 
point is a Hopf bifurcation. Furthermore, it is clear from 
Eq. (|16p that Oc depends only on the relative width pa- 
rameter /z, and it is also straightforward to establish that 
ttc increases with increasing /i, that is, a stronger coupling 
is necessary to overcome increasingly different values of 
71 and 72 (see lower panel of Fig. [S]). Note, however, 
that the dependence on /i implies that it is not just the 
difference in transition rates but the relative difference 
or percent difference relative to the mean transition rate 
that is the determining factor in how strong the coupling 
must be for synchronization to occur. A small-// expan- 
sion leads to an estimate of Oc to 0(/i^), 



12 + 3/i2 + \/3V(12 + Ai2)(4 + 3//2) 



(18) 



a result that exhibits these trends explicitly. The upper 
panel in Fig. [7] shows that this estimate is remarkably 
helpful even when fi^ is not so small. 

The frequency of oscillation of the synchronized system 
at the transition is given by w = lima_^a^ImA-|_. From 



Eq. ((T6)) it follows that u> depends on (71 +72) as well as 
fi. The small-/i expansion leads to the estimate 



LU = Im{X±)\a- 



7V3( 



7i+72)(4 + m'), (19) 



which works exceedingly well for all ^ (see Fig. [7]). 

To check the predictions of our linearization procedure, 
we numerically solve the nonlinear Af = 2 mean field 
equations. In agreement with the structure of the lin- 
earized eigenvalues, all components of P{t) synchronize 
to a common frequency as the phase boundary in (/i, a)- 
space) is crossed. Interestingly, the numerical solutions 
also give us insight into the amplitude of the oscillations; 
that is, they allow us to explore the relative "magnitude" 
of synchronization within the two populations. As we will 
see, the two populations indeed oscillate with the same 
frequency, but with amplitudes and "degrees of synchro- 
nization" that can be markedly different. Consider the 
order parameter r{t) given by 



N 



(20) 



Here (j>^ is the discrete phase 27r(fc-l)/3 for state k g 
{1, 2, 3} at site v. For phase transition studies, one would 
likely average this quantity over time in the long time 
limit, and also over independent trials. For our purposes 
here, though, we find the time-dependent form more con- 
venient. In the mean field case, where we solve for prob- 
abilities to be in each state, the order parameter is easily 



2.5r 




FIG. 8: (Color online) The components Pi, 71 and P2,7i (two 
lighter or brown curves), and Pi, 72 a-nd P2,72(two darker or 
blue curves), of the vector P{t) vs time for 7 = 1, A = 0.125, 
with a — 3.15, which is above the critical value ~ 3.02 pre- 
dicted by linearization. The left inset shows the order param- 
eter r{t) as it approaches its long-time limit. The right inset 
shows the frequency spectrum of a component of P{t). The 
spectrum has a dominant peak near ~ 4, and is expected 
to approach the frequency uj ~ 3.5 predicted by linearization 
as we approach the transition point a Uc- 
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FIG. 9: (Color online) The components Pi,-fi and ^2,71 (two 
lighter or brown curves), and -Pi, 72 and P2,72(two darker or 
blue curves), of the vector P{t) vs time for 7=1, A = 0.625, 
with a = 3.5, which is above the critical value Oc ~ 3.39 
predicted by linearization. The left inset shows the order 
parameter r{t) as it approaches its long-time limit. The right 
inset shows the frequency spectrum of a component of P{t). 
The spectrum has a dominant peak near lj ~ 4.4, and is 
expected to approach the frequency lu ~ 3.8 predicted by 
linearization as we approach the transition point a ^ ac- 
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FIG. 10: (Color online) The frequency uj of macroscopic os- 
cillations approaches the value predicted by linearization as 
a ttc. Here we have used 71 = 3.5 and 72 = 4.5 so that 
7 = 4 and A = 1. The critical coupling constant is flc = 3.06. 
From darkest or blue to lightest or brown: a = 3.07, 3.10, 
3.15, and 3.20. 



FIG. 11: (Color online) Time-averaged order parameter r in 
the long-time limit vs A for the individual oscillator popu- 
lations characterized respectively by the transition rate pa- 
rameter 71 (stars) and 72 (circles), and for the entire mixed 
array (squares). The insets show the time evolution of the 
probability vector components Pi, 71 and ^2,71 (lighter or 
brown curves) and Pi, 72 and ^2,72 (darker or blue curves) 
for widths 0.05 (upper inset), 0.5 (middle inset), and 0.9 
(lower inset). Some of the curves are not visible because they 
are so perfectly superimposed. While the degree of synchro- 
nization varies within each population, the critical width for 
de-syncrhonization is the same for both, as predicted by lin- 
earization. The coupling constant for all cases is a = 3.2 and 
the average transition rate parameter 7 = 1.5. 
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FIG. 12: (Color online) Same as Fig.[TT]but with 7 = 3.5 and 
widths 0.1 (upper inset), 1.3 (middle), and 1.8 (lower). 



calculated by writing the average in Eq. pp)) in terms of 
these probabilities rather than as a sum over sites. 

In subsequent figure captions we introduce the nota- 
tion 7 = (71 -1-72)72 (average transition rate parameter), 
and the difference A = I72 — 7i| (note that ^ ~ A/7). As 
shown in Figs.[8lfT2l the predictions of linearization accu- 
rately describe the onset of macroscopic synchronization 
and provide an estimation of the frequency of these oscil- 
lations near threshold (see Fig. fTU]) . Specifically, Figs. [5] 



and [S] show the macroscopic oscillations for coupling a 
above threshold. All the oscillators, regardless of their 
intrinsic transition rate parameter, oscillate exactly in 
phase, but the degree of synchronization is greater in the 
population with the larger (here 72), as evidenced by 
the unequal amplitude of the components of P{t) for the 
two populations. The "greater degree of synchroniza- 
tion" is also apparent in the order parameter r(t) shown 
in the insets, which is larger for the oscillators with the 



8 



higher intrinsic transition rate. These results support the 
notion that populations with higher transition rate pa- 
rameters in some sense synchronize more readily. The 
figures also show the frequency spectrum of any compo- 
nent of P{t). The peak occurs at the frequency of oscilla- 
tion of the synchronized array. As a — > Oc this frequency 
approaches the value ImA_|_ predicted by linearization, as 
shown in Fig. [TOl 

Figures [TT] and [T2] illustrate the sudden de- 
synchronization (at fixed a and 7) accompanying an in- 
crease in the difference A. This behavior is reminscent of 
that of the original Kuramoto oscillators, which become 
disordered as the width of the frequency distribution 
characterizing the population exceeds some critical value. 
The insets show the components of P{t) and confirm 
that both populations undergo the de-syncrhonization 
transition at the same critical value of the difference 
A. Comparing the two figures, we see that the sys- 
tem with a higher average transition rate parameter 
(Fig. [T^ can withstand a larger difference A before de- 
syncrhonization, again confirming our earlier observa- 
tions. 

One last point to consider is the relation between the 
frequency of oscillation of the synchronized array above 
ttc and the frequencies of oscillation of the two popula- 
tions if they were decoupled from one another. As cou- 
pling increases, the oscillation frequency uj moves closer 
to that of the population with the lower transition rate 
parameter. This is illustrated in Fig. [13] for the same 
parameters used in Fig. [TOl 

Finally, a visually helpful illustration of these behav- 
iors is obtained via a direct simulation of an array with 
a dichotomous population of oscillators. Since our os- 
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FIG. 14: (Color online) Long-time snapshots of a globally 
coupled system above and below threshold. In both cases, 
(71,72) = (0.5, 1.5). On the left, a — 3.5 < ac while on the 
right, a = 4.1 > Oc. In both cases, all units are globally 
coupled. For visualization purposes, the plot is arranged so 
that population 71 consists of the first 2500 units (the top), 
while population 72 consists of the second 2500 units (the 
bottom). Global synchrony emerges for a > ac- In addition, 
the population with the higher transition rate parameter is 
more synchronized. 



cillators are globally (all-to-all) coupled, the notion of a 
spatial distribution is moot, and for visulatization pur- 
poses we are free to arrange the populations in any way 
we wish. In Fig. [TH we display an equal number of 71 
and 72 oscillators and arrange the total polulation of 
N — 5000 so that the first 2500 have transition rate pa- 
rameter 71 and the remaining 2500 have transition rate 
parameter 72. In this simulation we have chosen 71 — 0.5 
and 72 = 1.5, so that 7 = A = 1. Although N = 5000 
is not infinite, it is large enough for this array to behave 
as predicted by our mean field theory. The left panel 
shows snapshots of the phases (each phase is indicated 
by a different color) when a < ac and the phases are ran- 
dom. The right panel shows the synchronized array when 
a > Qc- Clearly, all units are synchronized in the right 
panel, but the population with the higher transition rate 
parameter (lower half) shows a higher degree of synchro- 
nization (higher P{t)) as indicated by the intensity of the 
colors or the gray scale. 



B. A/" = 3 and J\f = A 



FIG. 13: (Color online) Frequency spectra of the numerical 
solutions to the mean field equations (A/" = 2) show that above 
critical coupling, a = 3.2 > Uc = 3.06, synchronization occurs 
at a frequency closer to the lower of the two population fre- 
quencies. Top left inset: Pi ,,j^=3.5 and Pi, 7=4. 5 when all units 
are globally coupled. Bottom right inset: the same curves for 
populations that are uncoupled from one another (but still 
globally coupled within each population). 



We can carry out this analysis, albeit not analytically 
(at least in practice), for any Af. We have explored the 
cases N' — 3 and 4. In both cases there appears to be only 
one pair of eigenvalues whose real parts can become pos- 
itive, suggesting that synchronization occurs all at once 
and not in one population at a time (Figs. 1151 and Il6p . 
This occurs no matter the distribution of the 3 or 4 tran- 
sition rate parameters. For example, in the A/" = 4 case 
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FIG. 15: Af — 3 case: The three pairs of complex conjugate 
eigenvalues are plotted in the complex plane for 2.9 < a < 4.5. 
The transition rate parameters are chosen to be (71, 72, 73) ~ 
(1, 2, 3). In the bottom panel the critical value of a is ac ~ 3.6. 



we have compared in some detail the cases where the four 
transition rates are equidistant and where they are pair- 
wise much closer than the separation between the highest 
and lowest. In both cases there is a single transition to 
synchrony, albeit not at exactly the same value of ac, 
indicating a more complex dependence on the transition 
rate parameter distribution than just via its mean and 
width. Furthermore, the basic trends of the dichotomous 
case broadly carry over, mainly in that the critical value 
Oc increases when the width of the distribution increases 
relative to the mean (as one would expect). On the other 
hand, the inclusion of more transition rates within a given 
range leads to a lowering of the critical coupling. Thus, 
for example, the mean transition rate 7 and the width 



A are the same in the cases shown in Figs. [5] and [15] 
(7 = A = 2), and yet Oc is higher in the former (3.95 for 
Af = 2) than in the latter (3.6 for Af = 3). Still, the mean 
and width of the distribution provide a rough qualitative 
assesmcnt of the behavior, particularly for the case of a 
uniform distribution, which we study below. 



V. UNIFORM DISTRIBUTION OF 
TRANSITION RATE PARAMETERS 

We now turn to globally coupled arrays where the tran- 
sition rate parameter g for each unit is chosen from a 
uniform distribution over a finite interval, (p{g). While 
it is difficult to make direct analytical progress in this 
general case, the earlier dimer analysis and the arrays of 
J\f = 2,3,4 different populations of units provide a frame- 
work for understanding the properties of these more gen- 
eral systems. In particular, the earlier results suggest 
that this "more disordered" system may also display a 
single transition to synchronization. To explore these 
and other features in more detail, we simulate N = 5000 
globally connected units characterized by the transition 
rate parameter distribution ip{g), and we make several 
observations. Firstly, wc do observe a single transition 
to macroscopic synchronization. Secondly, as suggested 
by the dichotomous case, synchronization appears more 
readily (that is, for a lower value of a) if the distribu- 
tion (/)(g) has a larger mean and smaller width. When 
the mean and width are varied independently, the qual- 
itative trends from the dichotomous case are observed 
here as well. Thirdly, while synchronization in this sys- 
tem is again governed primarily by the mean and width 
of the distribution (p{g), the critical value Oc is consid- 
erably lower than that of the finite J\f systems with the 
same mean and width (as expected). 

Two examples of our simulation results are shown in 
Figs. [T71 and [T51 In Fig.[T7]we present the first two com- 
ponents of the 3-dimensional vector P{t) whose compo- 
nents Pi (t) represent the probability that all units of the 
entire synchronized array are in state i. The probabili- 
ties Pi (t) and P2 (t) oscillate in time with essentially con- 
stant amplitude and a constant relative phase, indicat- 
ing global synchronization. The upper left inset shows 
the order parameter r{t) and the upper right inset the 
time resolved snapshot of the system, both indicating a 
high degree of synchronization. Note that the coupling 
parameter a = 3 in the figure is below the critical value 
flc = 3.2 for the dichotomous case with the same mean 
and width. 

Figure [T8l shows the steady state time-averaged order 
parameter r at constant a as the width of the (f>{g) distri- 
bution is increased for a fixed mean. Similar to the J\f = 2 
population case, synchronzation is destroyed as the width 
eclipses some critical value, and that value increases as 
the mean of the distrubtion increases. In Fig. [11] we plot 
the data from Fig [TS] as a function of the relative width 
parameter Recalling that for the dichotomous array as 
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FIG. 16: A/" = 4 case: The four pairs of complex conju- 
gate eigenvalues are plotted in the complex plane for 2.9 < 
a < 4.5. The transition rate parameters are chosen to be 
(71,72,73,74) = (1,2,3,4). In the bottom panel, the critical 
value of a is ac ^ 3.75. 



FIG. 17: (Color online) The probability that the synchronized 
array is in state 1 (lighter or brown) and state 2 (darker or 
blue) as a function of time for a uniform distribution ip{g) 
on the interval [1.5, 2.5] and coupling parameter a = 3. In- 
sets show the order parameter r{t) as well as time resolved 
snapshots of the system. 
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FIG. 18: (Color online) As the width of the (j>{g) distribution 
increases, a critical width is reached beyond which synchro- 
nization is destroyed. The coupling is chosen to be a = 3.2, 
and the four curves represent the steady state, time-averaged 
order parameter for distributions with different means. As 
the mean of the 4>{g) distribution increases, the transition to 
disorder occurs at a greater width. The insets at the right 
show the long-time behavior of an entire population of mean 
transition rate parameter 3.5 (corresponding to the triangle 
order parameter data) and widths of 0.6, 4.0, and 6.2. 



well as for the dimer synchronzation at a given a depends 
only on fi, we might expect that the transition point fic 
(at constant a) is not significantly mean-dependent, even 
when there is a distribution of transition rate parame- 
ters. In fact, we can see that the curves approximately 
collapse onto one curve, indicating that the relative width 
fi provides a useful control parameter for predicting syn- 
chronization. Hence, the predictions of the linearization 
analysis for the Af = 2 case provides qualitatively insight 
into the behavior of the disordered population. 
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FIG. 19: The data of Fig. 1181 against the relative width fi. 



VI. DISCUSSION 

We have presented a discrete model for globally cou- 
pled stochastic nonlinear oscillators with a distribution of 
transition rate parameters. Our model exhibits a range 
of interesting dynamical behavior, much of which mimics 
the qualitative features of the canonical Kuramoto os- 
cillator but with a mathematically and numerically 
considerably more tractable model. Since our phase vari- 
able is discrete (whereas the phase variable in the canon- 
ical problem is continuous), a distribution of Af different 
transition rates in our array leads to a set of 2Af coupled 
nonlinear ordinary differential equations instead of a sin- 
gle partial differential equation for the probability distri- 
butions of interest. Linearization of our model around 
the critical point leads to a problem which at least for 
small M (specifically, for the dichotomous disorder case) 
becomes analytically tractable. Distributions involving a 
large finite number of transition rate parameters, while 
not easily amenable to analytic manipulation even upon 
linearization, reduce to a simple matrix algebra problem. 
For any distribution of transition rate parameters, even 
continuous, the model is in any case readily amenable to 
numerical simulation. 

Our most salient conclusion is that such disordered 



globally coupled arrays of oscillators, even in the face of 
transition rate parameter disorder, undergo a single tran- 
sition to macroscopic synchronization. Furthermore, we 
have shown that the critical coupling ac for synchroniza- 
tion depends strongly (but not exclusively) on the width 
A and mean 7 of the transition rate parameter distribu- 
tion, specifically via the relative width fi = A/7. This 
general feature is already apparent in the synchronization 
behavior of a dimer of two oscillators with transition rate 
parameters 71 and 72. An infinite array of two popula- 
tions of oscillators, one with transition rate parameter 71 
and the other with 72 , displays a Hopf bifurcation, with 
flc determined solely by /x. While a quantitative predic- 
tion of synchronization on the basis of the relative width 
is not possible in all cases, it does determine qualitative 
aspects of the transition for more complex transition rate 
parameter distributions. We have explored this assertion 
for arrays with Af = 2, 3, and 4 and with a uniform 
distribution of transition rates over a finite interval, and 
expect it be appropriate for other smooth distributions 
as well. 

A number of further avenues of investigation based on 
our stochastic three-state phase-coupled oscillator model 
are possible. For example, we could explore the effects of 
transition rate disorder in locally coupled arrays whose 
behavior we have fully characterized for identical oscil- 
lators [1, . It would be interesting to explore the con- 
sequences of disorder in the coupling parameter a. Fi- 
nally, we note that a two-state version of this model 
(which of course does not lead to phase synchroniza- 
tion as discussed here) has recently been shown to accu- 
rately capture the unique statistics of blinking quantum 
dots [l3| . Such wider applicability of the model, together 
with its analytic and numerical tractability, clearly opens 
the door to a number of new directions of investigation. 
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